# Load Necessary Packages
library(APCI)
library(ggplot2)
library(tidyr)
library(dplyr)
library(ggpubr)

# Load Datasets
# setwd("")
data <- read.csv("IE data.csv",header=T,sep=",")

# Variable names
data$monarchism <- data$mo.rchism

# Figure 1
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
data.overtime.1 <- data %>% group_by(period) %>% summarise(m=weighted.mean(monarchism,weight,na.rm=T),a=weighted.mean(abolition,weight,na.rm=T))
data.overtime <- data.frame(period=rep(data.overtime.1$period,2),value=c(data.overtime.1$m,data.overtime.1$a),series=factor(c(rep("Monarchism",21),rep("Abolition",21)),levels=c("Monarchism","Abolition")))
p.overtime <- ggplot(data.overtime,aes(x=period,y=value,group=series,color=series)) + geom_point() + geom_line() + theme_classic() + scale_color_manual(name="Variable",values=cbbPalette) + ylim(0,1) + ylab("Support") + xlab("Year")
# ggsave("Figure 1.png",p.overtime,device="png",width=6,height=4,units="in")

# Run Baseline APCs
## Monarchism
mod.baseline.m <- APCI::temp_model(data=data,outcome='monarchism',age='agegroup',period='periodgroup',covariate=c('country'),family='gaussian')
cohort.baseline.m <- APCI::cohortdeviation(mod.baseline.m$A,mod.baseline.m$P,mod.baseline.m$C,model=mod.baseline.m$model,covariate=c('country'))
## Abolition
mod.baseline.a <- APCI::temp_model(data=data,outcome='abolition',age='agegroup',period='periodgroup',covariate=c('country'),family='gaussian')
cohort.baseline.a <- APCI::cohortdeviation(mod.baseline.a$A,mod.baseline.a$P,mod.baseline.a$C,model=mod.baseline.a$model,covariate=c('country'))

# Figure 2
## Grand Means
weighted.mean(data$monarchism,w=data$weight,na.rm=T) # Monarchism=0.70
weighted.mean(data$abolition,w=data$weight,na.rm=T) # Abolition=0.09
## Monarchism
f2data.m <- data.frame(age=c(15,20,25,30,35,40,45,50,55,60,65,70,75,80,85),
                       effect=summary(mod.baseline.m$model)$coefficients[4:18,1],
                       se=summary(mod.baseline.m$model)$coefficients[4:18,2])
f2data.m$lci <- f2data.m$effect-(1.96*f2data.m$se)
f2data.m$uci <- f2data.m$effect+(1.96*f2data.m$se)
f2data.m$sig <- factor(ifelse(f2data.m$lci<0 & f2data.m$uci>0,"No","Yes"))
## Abolition
f2data.a <- data.frame(age=c(15,20,25,30,35,40,45,50,55,60,65,70,75,80,85),
                       effect=summary(mod.baseline.a$model)$coefficients[4:18,1],
                       se=summary(mod.baseline.a$model)$coefficients[4:18,2])
f2data.a$lci <- f2data.a$effect-(1.96*f2data.a$se)
f2data.a$uci <- f2data.a$effect+(1.96*f2data.a$se)
f2data.a$sig <- factor(ifelse(f2data.a$lci<0 & f2data.a$uci>0,"No","Yes"))
## Joint Figure
f2.m <- ggplot(f2data.m,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Age, Monarchism") + theme(plot.title = element_text(hjust = 0.5))
f2.a <- ggplot(f2data.a,aes(x=age,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Age") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Age, Abolition") + theme(plot.title = element_text(hjust = 0.5))
f2 <- ggpubr::ggarrange(f2.m,f2.a,nrow=1)
# ggsave("Figure 2.png",f2,device="png",width=6,height=4,units="in")

# Figure 3
## Monarchism
f3data.m <- data.frame(cohort=seq(1890,2000,by=5),
                       effect=as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average)),
                       lci=as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average_se))),
                       uci=as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.baseline.m$cohort_average$cohort_average_se))))
f3data.m$sig <- ifelse(f3data.m$lci<0 & f3data.m$uci>0,"No","Yes")
## Abolition
f3data.a <- data.frame(cohort=seq(1890,2000,by=5),
                       effect=as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average)),
                       lci=as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average))-(1.96*as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average_se))),
                       uci=as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average))+(1.96*as.numeric(as.character(cohort.baseline.a$cohort_average$cohort_average_se))))
f3data.a$sig <- ifelse(f3data.a$lci<0 & f3data.a$uci>0,"No","Yes")
## Joint Figure
f3.m <- ggplot(f3data.m,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Average, Monarchism") + theme(plot.title = element_text(hjust = 0.5))
f3.a <- ggplot(f3data.a,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Average, Abolition") + theme(plot.title = element_text(hjust = 0.5))
f3 <- ggpubr::ggarrange(f3.m,f3.a,nrow=1)
# ggsave("Figure 3.png",f3,device="png",width=6,height=4,units="in")

# Figure 4
## Monarchism
f4data.m <- data.frame(cohort=seq(1885,1995,by=5),
                       effect=as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope)),
                       lci=as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope_se))),
                       uci=as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.baseline.m$cohort_slope$cohort_slope_se))))
f4data.m$sig <- ifelse(f4data.m$lci<0 & f4data.m$uci>0,"No","Yes")
## Abolition
f4data.a <- data.frame(cohort=seq(1885,1995,by=5),
                       effect=as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope)),
                       lci=as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope))-(1.96*as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope_se))),
                       uci=as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope))+(1.96*as.numeric(as.character(cohort.baseline.a$cohort_slope$cohort_slope_se))))
f4data.a$sig <- ifelse(f4data.a$lci<0 & f4data.a$uci>0,"No","Yes")
## Joint Figure
f4.m <- ggplot(f4data.m,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Slopes, Monarchism") + theme(plot.title = element_text(hjust = 0.5))
f4.a <- ggplot(f4data.a,aes(x=cohort,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Cohort") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Cohort Slopes, Abolition") + theme(plot.title = element_text(hjust = 0.5))
f4 <- ggpubr::ggarrange(f4.m,f4.a,nrow=1)
# ggsave("Figure 4.png",f4,device="png",width=6,height=4,units="in")

# Figure A1
## Monarchism
fa1data.m <- data.frame(period=c(1990,1995,2000,2005,2010,2015,2020),
                        effect=summary(mod.baseline.m$model)$coefficients[19:25,1],
                        se=summary(mod.baseline.m$model)$coefficients[19:25,2])
fa1data.m$lci <- fa1data.m$effect-(1.96*fa1data.m$se)
fa1data.m$uci <- fa1data.m$effect+(1.96*fa1data.m$se)
fa1data.m$sig <- factor(ifelse(fa1data.m$lci<0 & fa1data.m$uci>0,"No","Yes"))
## Abolition
fa1data.a <- data.frame(period=c(1990,1995,2000,2005,2010,2015,2020),
                        effect=summary(mod.baseline.a$model)$coefficients[19:25,1],
                        se=summary(mod.baseline.a$model)$coefficients[19:25,2])
fa1data.a$lci <- fa1data.a$effect-(1.96*fa1data.a$se)
fa1data.a$uci <- fa1data.a$effect+(1.96*fa1data.a$se)
fa1data.a$sig <- factor(ifelse(fa1data.a$lci<0 & fa1data.a$uci>0,"No","Yes"))
## Joint Figure
fa1.m <- ggplot(fa1data.m,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Period, Monarchism") + theme(plot.title = element_text(hjust = 0.5))
fa1.a <- ggplot(fa1data.a,aes(x=period,y=effect,color=sig)) + geom_point() + theme_classic() + geom_errorbar(aes(ymin=lci,ymax=uci),lwd=0.5,width=0) + geom_hline(yintercept=0,linetype="dashed",color="red") +xlab("Period") + ylab("Difference from Grand Mean") + scale_color_manual(values=c("Grey50","Black"),guide="none") + ggtitle("Period, Abolition") + theme(plot.title = element_text(hjust = 0.5))
fa1 <- ggpubr::ggarrange(fa1.m,fa1.a,nrow=1)
# ggsave("Figure A1.png",fa1,device="png",width=6,height=4,units="in")
